#include "include/Element.hpp"
#include <iostream>
#include <stdlib.h>
#include <cmath>
#include <array>

extern int degree_of_basis_function;
extern int element_shape;

Element::Element(unsigned long idx, std::initializer_list<Point> points, bool boundary)
{
    index = idx;
    for(auto iter=points.begin(); iter < points.end(); iter++)
    {
        vertex.push_back(*iter);
    }
    marker = boundary;
}

Point Element::operator[](unsigned int idx)
{
    return vertex[idx];
}

std::ostream& operator<<(std::ostream& os, Element e)
{
    os << "Element "<< e.index << " has points: " << std::endl;
    for(auto iter=e.vertex.begin(); iter<e.vertex.end(); iter++)
    {
        os << (*iter);
    }
    if(is_boundary(e) == 1)
    {
        os << "on the boundary" << std::endl;
    }
    else
    {
        os << "inside the domain" << std::endl;
    }
    return os;
}

unsigned long id(const Element& e)
{
    return e.index;
}

bool is_boundary(const Element& e)
{
    return e.marker;
}

double Jacobi(const Element& e)
{
    double rst = 0.0;
/*assemble matrix
 *  (b2  b3)
 *  (c2  c3)
 */
    double b2 = 0.0;
    double b3 = 0.0;
    double c2 = 0.0;
    double c3 = 0.0;
    if(element_shape == 3)
    {
        Point p1 = e.vertex[0];
        Point p2 = e.vertex[1];
        Point p3 = e.vertex[2];
        b2 = p3[1] - p1[1];
        b3 = p1[1] - p2[1];
        c2 = p1[0] - p3[0];
        c3 = p2[0] - p1[0];
        rst = b2 * c3 - b3 * c2;
    }
    else
    {
        std::cerr << "Qualateral grid is not supported yet!" << std::endl; 
    }
    if(rst == 0)
    {
        std::cerr << "Invalid element because area is zero." << std::endl;
        abort();
    }
    return rst;
}

// double Element::integral(int i, int j,  Region& r)
// {
    // double J = Jacobi(*this);
    // double rst = 0.0;
    // for(auto iter=r.nodes.begin(); iter < r.nodes.end(); iter++)
    // {
        // rst += (*iter).get_weight() * r.basis_gradient[i][j]( (*iter)[0], (*iter)[1] );
    // }
    // rst = rst * J;
    // return rst;
// }

double Element::integral(std::function<double(double, double)> f, Region& r)
{
    double rst = 0.0; 
    for(auto n=r.nodes.begin(); n<r.nodes.end(); n++)
    {
        std::array<double, 2> coord = l2g((*n)[0], (*n)[1]);
        double x = coord[0];
        double y = coord[1];
        rst += f(x,y) * ((*n).get_weight());
//        std::cerr << "global integral point(" << x << ", " << y << ") with weight=" << ((*n).get_weight()) << " value=" << temp << std::endl;
    }
//    std::cout << "final integral evaluate=" << rst << std::endl;
    rst = rst * abs(Jacobi(*this)) / 2.0;
    return rst;
}

std::array<double, 2> Element::l2g(double xi, double eta)
{
    double x = vertex[0][0] * (1.0 - xi - eta) + vertex[1][0]* xi + vertex[2][0] * eta;
    double y = vertex[0][1] * (1.0 - xi - eta) + vertex[1][1]* xi + vertex[2][1] * eta;
    std::array<double, 2> rst = {x, y};
    return rst;
}

std::array<double, 2> Element::g2l(double x, double y)
{
    std::array<double, 2> rst = {0.0, 0.0};
    if(element_shape == 3)
    {
        double x1 = vertex[0][0];
        double x2 = vertex[1][0];
        double x3 = vertex[2][0];
        double y1 = vertex[0][1];
        double y2 = vertex[1][1];
        double y3 = vertex[2][1];
        
        rst[0] = ((x-x1)*(y3-y1) - (y-y1)*(x3-x1)) / ((x2-x1)*(y3-y1) - (y2-y1)*(x3-x1));
        rst[1] = ((x-x1)*(y2-y1) - (y-y1)*(x2-x1)) / ((x3-x1)*(y2-y1) - (y3-y1)*(x2-x1));
    }
    else
    {
        std::cerr << "only triangle elements with P1 is feasible" << std::endl;
    }
    return rst;
}

void Element::add_point(Point& np)
{
    vertex.push_back(np);
}
